Fractal fronts in fractal fractures: large and small-scale structure 
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Abstract 

The evolution and spatial structure of displacement fronts in fractures with self-affine rough walls 
are studied by numerical simulations. The fractures are open and the two faces are identical but 
shifted along their mean plane, either parallel or perpendicular to the flow. An initially flat front 
advected by the flow is progressively distorted into a self-affine front with Hurst exponent equal to 
that of the fracture walls. The lower cutoff of the self-affine regime depends on the aperture and 
lateral shift, while the upper cutoff grows linearly with the width of the front. 

PACS numbers: 47.55Mh, 05.40.-a, 92.40.Kf, 47.53. +n 



1 



Fractures are ubiquitous in the Earth's upper crust and in many materials of engineer- 
ing interest, and their size varies over a broad range of length scales, from several meters 
in the case of faults down to the sub-micron scale of micro-cracks. The geometry of the 
fracture obviously influences any transport process within it, and in particular the hydraulic 
properties of rock fractures, crucial in such applications as hydrocarbon recovery, subsurface 
hydrology, and earthquake prediction. In these situations, transport generally involves flow 
through entire networks of fractures, and often through a porous matrix as well, but fre- 
quently the flow is dominated by channeling through particular fractures (jj,!^, and therefore 
there is considerable interest in fully understanding flow through an individual fracture. 
A key feature characterizing the pore space geometry of fractured materials is the exper- 
imental observation that the roughness of fracture surfaces is statistically described as a 
self-affine fractal structure over an appreciable range of length scales P, |^ . This observation 
holds for both man-made and natural fractures [8|, |9|, and for a broad variety of 

materials and length scales. We consider rock fracture surfaces without overhangs, whose 
height can be described through a single valued function z{x,y), where {x,y) are Cartesian 
coordinates in the mean plane of the surface. Self-affinity means that the surface remains 
statistically invariant under the scale transformation z{\x, Xy) = X''z{x,y), where the Hurst 
exponent ( characterizes the roughness of the surface, in that the fluctuation of the surface 
heights over a length L is given by cr^^L) = i{L/tj^. Here ^ is the topothesy, defined as 
the horizontal distance over which fiuctuations in height have a rms slope of one [9]. Ex- 
perimental values of C are found to be close to 0.8, for most of the materials considered and 
indepe„de. of the .actu.e .ode g B B Q. w.h the exception of .ate.a.s d,sp.ay- 
ing intergranular fractures, such as some sandstone rocks, for which it has been found to 
be close to 0.5 [10]. The same exponent C, describes the scaling behavior of the two-point 
height correlation function, and therefore, one expects the fractal geometry of the surfaces 
to induce spatially correlated velocity fiuctuations. 

A particularly sensitive measurement of the correlated character of the velocity fiuctua- 
tions is the shape of an advected front of tracer particles, because each particle position 
is an integral over the velocity field along its path. Recent displacement experiments have 
found a striking approximate correlation between tracer motion and surface geometry - the 
shape of the front was found to be self-affine, with an exponent close to the Hurst exponent 
C = 0.8 of the fracture wall The objective of this paper is to establish this connection 
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FIG. 1: Fronts obtained when the lateral to the flow (top; dx = 0.4 mm) and 

parallel to the flow (bottom; dy = 0.4 mm). The black lines correspond to a fracture with Hurst 
exponent = 0.5 and the gray lines correspond to ^ = 0.8. The inserts display 4 mm x 4 mm 
sections of the fronts. 



in quantitative detail, and determine its range of validity and further implications, using 
Lattice-Boltzmann (LB) numerical simulations We consider the motion of passive trac- 
ers in single fractures, which are modeled as the opening between two perfectly matching 
self-affine rough surfaces. We assume that there is no contact between the two fracture 
surfaces, and that the walls are chemically inert, completely rigid, and non-porous. As in 
the experiments motivating this work, we focus on the case of small Reynolds numbers for 
which inertia effects can be neglected, and Peclet number sufficiently high to ignore molec- 
ular diffusion. 

A self-affine surface is first generated numerically, using a Fourier synthesis method 



and the fracture pore space is the region between one such surface and a suitably shifted 
rephca. First, the two surfaces are separated by a fixed distance h, in the direction normal 
to the mean plane of the fracture, and then a lateral shift of the upper surface is introduced. 
The effects of the latter lateral shift, meant to reflect subsurface motion during or subse- 
quent to fracture, are an important practical issue. The maximum allowed shift corresponds 
to the occurrence of the first contact point between the two surfaces, so that the fracture is 
open over its full area. Two orthogonal shift directions along the mean plane of the fracture 
are considered; dy, in the direction of the mean flow, and dx, perpendicular to it. The free 
space between the two surfaces is occupied by a cubic lattice of fixed lattice spacing, with 
1024 X 1024 X 20 sites, which is large enough to accurately capture the three-dimensionality 
of the velocity field |l6|. 

In order to generate numerically fractures that are comparable to those used in the ex- 
periments, we need to match the two relevant length scales, the topothesy and the mean 
aperture, as well as the Hurst exponent. Thus, the mean aperture is set to 1 mm, corre- 
spondingto a fracture size L = 51.2 mm {6 = 0.05 mm), equal to one fourth of the fracture 
used in \L^. In order to have the same surface roughness, the amplitude of average fluctu- 
ations in surface height over a length of 1 mm is set to 0.25 mm. Two different values of 
the Hurst exponent ( are used here: 0.8 and 0.5. The first value corresponds to the mea- 



sured Hurst exponent of the material used in 



12j , and the second one is typical of fractured 



materials with intergranular effects. Considerable height variations occur over a distance 
comparable to the aperture, az{h) = 0.25 mm and, therefore, the lubrication approximation, 
which assumes a Poiseuille velocity profile across the aperture, fails to describe the fiow field 
inside the fractures 111. Therefore, the three-dimensional fiow field in the pore space of 

n 

the fractures is calculated, for different values of the lateral shift, using a LB method |17| . 
Periodic boundary conditions are used in the x and y directions, including the inflow and 
outflow, and the Reynolds number is small (~ 0.2). 

n 

In the displacement experiments presented in Ref. a steady radial flow of a clear 

fluid is flrst established in a transparent fracture model, cast in epoxy from natural granite. 
The same fluid, with a small concentration of a dye, is then injected and the spreading front 
of the dye is recorded. Several mechanisms contribute to the dispersion of tracer particles, 
among them molecular diffusion and Taylor dispersion. However, the most relevant process 
that drives the spreading within fractures is the velocity fluctuations induced by the sur- 
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TABLE I: Statistical properties of the front in mm displayed in Fig. ^ (y) is the mean position of 
the tracers in the flow direction, a = {{y — (y))'^)^^'^, and e = ymax — Vmin 
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FIG. 2: Mean velocity of the tracers normalized by the velocity for zero shift as function of the 
shift, o and o are for flow along the shift direction within fractures of wall roughness C, = 0.5 and 
0.8. □ and A correspond to flow normal to the shift in fractures with Q = 0.5 and 0.8. 

face roughness |12j]. This mechanism, known as geometric or macro-dispersion, is linear in 
the velocity, and therefore leads to an external front that is independent of the flow rate 
for a flxed injected volume. The front is thus essentially determined by the two dimen- 
sional geometrical disorder of the velocity fleld in the mean fracture plane. Comparisons 
between experimental and numerical dispersion results can thus be achieved by considering 
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the spreading of tracer particles moving in the flow field resulting from averaging the three 
dimensional velocities over the gap of the fracture. 

Initially in the computations, a tracer particle is located at each of the 1024 lattice sites 
corresponding to the injection side of the fracture. Each particle then propagates in the gap 
averaged velocity field and the velocity fluctuations due to the effective aperture variation 
will distort the initially fiat front. Figure ^ shows typical invasion fronts obtained after 
the tracer traverses roughly 80% of the system. An obvious qualitative feature is that the 
front is rougher when the relative shift of the fracture surfaces is normal to the mean flow. 
Quantitatively, these small-scale features vary significantly with both the Hurst exponent 
and the shift direction, and in particular, the front contains more small scale details when 
the fracture wall roughness is decreased. The basic numerical properties of these fronts are 
entirely consistent with their visual appearance, as shown in Table HI We will show presently 
that these fronts are in fact fractals, over a range of length scales. In addition, in Fig. |21 
we show that there is a remarkable anisotropy in the average propagation velocity of the 
front, in that a relative shift perpendicular to the mean flow has a negligible effect, whereas 
a shift along the flow produces a systematic decrease in the permeability. Qualitatively, this 
behavior may be understood from the observation that transverse shifts simply replace one 
preferred flow channel with another of comparable hydraulic resistance while longitudinal 
shifts tend to block the flow, but the systematic variation, or lack of it, seen here is quite 
surprising [3]. 

The experimental results presented in were analyzed for their fractal dimension D 
using the average mass method j^, in which a constant density of tracer is assumed for the 
invasion front and the average number of tracer particles, (M(r)), located at a distance less 
than r from a random origin on the contour is computed. For a fractal front, a power law 
relationship (M(r)) oc is expected, where D is the fractal dimension ranging between 1 
and 2, and if the front is furthermore assumed to be self-affine, one would identify (f = 2 — D. 
The result found there was Cf = 0.74 ± 0.03, close to the Hurst exponent of the fractures, 

n 

( = 0.8. If we apply the same analyzis to the numerical results here, we obtain [1^ Q = 
0.82 ± 0.01 when ( = 0.8 and (j = 0.5 ± 0.15 for ( = 0.5, for length scales between certain 
upper and lower cutoffs. The lower cutoff is of the order the mean aperture, and independent 
of the direction and magnitude of the lateral shift, as were the dimensions (f themselves. 
However, the upper cutoff is shift- dependent, and difficult to interpret. An additional and 



6 




1 
Iog^o(a5) 



FIG. 3: Variation of logio{{\W[y] (a, as function of logiQ{a6). Circles and triangles correspond 

to C = 0.8 and ^ = 0.5 respectivley. Solid and open symbols correspond to dy = 0.4 and dx = 0.4 
respectively. The two data sets are shifted vertically for convenience. The solid line has a slope of 
2.6 while the dashed line has a slope of 2. 

severe shortcoming of this procedure is that it does not directly test self-affinity. 

A more appropriate method for self-affine structures is the average wavelet coefficient 
(AWC) analysis where a one dimensional front y(x) is transformed into the wavelet 
domain as 

Wiy]{a,b) = 4= / '^lbix)yi^)d^, (1) 



where ipa,b{x) is obtained from the analyzing wavelet tp, in our case a Daubechies wavelet 
jigj ]. via rescaling and translation, ipa,b{x) = ip{{x — b)/a). The AWC measures the average 
"energy" present in the front at a given scale, defined as the arithmetic average of |iy[y](a, 6)p 
over all possible locations b, and for a statistically self-affine profile with exponent C/, it scales 
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FIG. 4: /o5fio((|VF[j^] (a, as function of the normalized scale parameter Iogio{a6/(T). The data 

set are shifted vertically to make comparisons easier. They correspond to fronts obtained after 
different mean traveled distances (y), in a fracture of roughness ^ = 0.8 and with walls shifted by 
dx = 0.4mm. The insert shows the variation of o" as a function of the mean traveled distance. 

as {\W[y]{a,b)\'^)b ~ a^^^+^. Typical results of applying the AWC method to the numerical 
fronts shown in Fig. ^are given in Fig. El and we see a robust scaling behavior, depending 
only on the roughness exponent of the fracture surfaces and independent of the relative shift 
and the invasion time. We find (f = 0.82 ± 0.1 for surfaces with roughness ( = 0.8 and 
(f = 0.5 ±0.15 for surfaces with ( = 0.5, consistent with the values found using the average- 
mass method. These results are in good agreement with previous experimental studies of 
tracer dispersion and with a theoretical analysis based on a perturbative approach [2^ . 

As usual, one sees scaling behavior only between certain cutoff length scales. The origin 
of the upper cutoff can be understood by considering the "energy" as a function of a6/a, 
where a is the width of the front. In Fig. El data for the case ( = 0.8 and dx = 0.4 are 
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FIG. 5: Variation of logio{{\W[y](a,b)\'^)h) as function of logio{a6/a). The data sets are shifted 
vertically to make comparisons easier. They correspond to fronts obtained for the same average 
traveled distance {y) = 40.7 mm and for the same fracture ((" = 0.8), but for different values of the 
lateral shift. The insert displays the corresponding fronts for increasing shifts from top to bottom. 

plotted at different times, showing that the upper cutoff for scahng behavior is proportional 
to the width of the front. The width itself grows roughly linearly as a function of time, or 
equivalently with the mean distance traversed by the front, as seen in the inset of Fig. El so 
the upper limit of self-affine behavior increases with time. To understand the lower cutoff, 
it is instructive to compare fronts at the same mean traversed distance for different values 
of the shift. In Fig. Elwe display data for the ( = 0.8 fracture with different values of lateral 
shift, at the same transit time; the lower cutoff is seen to increase with increasing shift. The 
connection between scaling breakdown and shift can be elucidated by simply examining the 
front shapes in these four cases, as shown in the inset to Fig. El An increasing shift tends 
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to produce stagnant low-velocity regions where the front halts, and these regions destroy 
the self-similarity. The AWC method is particularly well-suited to address the lower cutoff 
issue, because the stagnant tails of the front contribute to the wavelet "energy" strongly at 
low values of a. In the average-mass method, in contrast, they contribute predominantly at 
large r whereas the lower cutoff is on the order of the mean aperture, independent of the 



surface exponent and any lateral shift 
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There is a an obvious intuitive appeal to the idea that a fractal geometry will somehow 
impose itself on a suitably flexible object propagating through it. In this paper we have 
used numerical simulations to expand upon the experimental observation that a (high Peclet 
number) displacement front in a self-affine fracture behaves in this way. We have presented 
compelling evidence that self-affine fronts indeed appear in this situation, for a wide range 
of parameters, and have further determined the relevant range of length scales for scaling 
behavior. This work raises a number of further intriguing questions, including (1) The 
behavior of higher moments of the front geometry; (2) The interplay between the self-affine 
character and the strong anisotropic behavior in ra dm/ injection flows; (3) The effects of finite 
Peclet numbers, which would presumably alter the scaling cutoffs lengths, on hydrodynamic 
dispersion. We hope to address these and other questions in future work. 
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